Objetivo: explorar la importancia de evaluar productos de precipitación, consideraciones clave, funciones objetivo utilizadas y el preprocesamiento de datos necesario antes de la evaluación.
Los productos de precipitación son fundamentales para comprender y monitorear el ciclo del agua y su impacto en diversos sectores como la agricultura, hidrología, gestión de recursos hídricos, prevención de desastres, entre otros.
La evaluación de la calidad y precisión de estos productos es crucial para garantizar la confiabilidad de los datos y su utilidad en la toma de decisiones y aplicaciones relacionadas.
Resolución espacial y temporal (Benedict et al., 2019)
El desempeño de los productos (Baez-Villanueva et al., 2020)
Hora de reporte de las estaciones (Beck et al., 2019)
Las funciones objetivo son métricas utilizadas para evaluar el desempeño de los productos de precipitación. Algunas funciones objetivo comunes son:
Adicionalmente, queremos analizar el desempeño de los productos en detectar ciertas intensidades de precipitación, para esto podemos usar funciones objetivo categoricas.
Para una descripción mas detallada de estas funciones objetivo categóricas puedes leer Zambrano-Bigiarini et al., (2017) y Baez-Villanueva et al., (2018)
En este módulo, utilizaremos el índice KGE’ (Kling-Gupta Efficiency) ya que es capaz de descomponer el rendimiento total en tres componentes distintos:
\[ \textrm{KGE'} = 1 - \sqrt{ (r-1)^2 + (\beta-1)^2 + (\gamma-1)^2 } \]
\[ r = \frac{ \sum_{i=1}^n{(O_i - \bar{O}) (S_i - \bar{S}) } }{\sqrt{\sum_{i=1}^n{(O_i - \bar{O})^2}}\sqrt{\sum_{i=1}^n{(S_i - \bar{S})^2}}} \]
\[ \beta = \frac{\mu_{s}}{\mu_{o}} \]
\[ \gamma = \frac{CV_{s}}{CV_{o}} = \frac{\sigma_{s}/\mu_{s}}{\sigma_{o}/\mu_{o}} \]
Antes de la evaluación, es necesario realizar ciertos pasos de preprocesamiento en los datos de precipitación.
Para este ejercicio, utilizaremos datos diarios de estaciones de la base de datos Global Historical Climatology Network daily (GHCNd) para el sur de África. La evaluación se realizará durante el periodo 1990-1992 y, para esto, también utilizaremos datos diarios de MSWEPv2.8 y CHIRPSv2.
Paso 1: Primero cargaremos las series temporales observadas y los metadatos. Para esto, almacenemos las rutas a los archivos:
observations_path <- 'C:/Users/Data/timeseries_1990-1992.zoo'
metadata_path <- 'C:/Users/Data/metadata.csv'Después, utilizaremos la función read.zoo cargando el paquete hydroGOF para leer las series temporales, y la función read.csv para leer los metadatos.
observations <- read.zoo(observations_path, header = TRUE)
metadata <- read.csv(metadata_path)
dim(observations)
## [1] 1096 1377
dim(metadata)
## [1] 1377 6Podemos ver que tenemos un total de 1377 estaciones. Pero, ¿cuántos datos faltantes tienen estas estaciones?
observations[1:5, 1:5]
## BC000068026 BC000068032 BC000068226 BC000068328 BC004612080
## 1990-01-01 0.0 0 0 0 0
## 1990-01-02 1.0 0 0 0 0
## 1990-01-03 0.5 0 0 0 0
## 1990-01-04 34.0 0 0 0 0
## 1990-01-05 14.5 0 0 0 0
head(metadata)
## ID LATITUDE LONGITUDE ELEMENT YEAR_START YEAR_END
## 1 BC000068026 -18.367 21.850 PRCP 1932 1992
## 2 BC000068032 -19.983 23.417 PRCP 1921 2022
## 3 BC000068226 -24.017 21.883 PRCP 1958 1992
## 4 BC000068328 -26.050 22.450 PRCP 1939 1992
## 5 BC004612080 -26.467 20.617 PRCP 1960 2022
## 6 BC005452830 -25.220 25.670 PRCP 1921 1989Paso 2: Analizar el número de datos faltantes en cada estación. Para esto, vamos a crear una función básica:
Ahora la aplicaremos a cada una de las series temporales (columnas del objeto zoo).
Podemos generar un histograma para visualizar la distribución de los datos faltantes:
Seleccionemos un umbral de 0.2 como límite de datos faltantes. Es decir, conservaremos las estaciones que tengan menos del 20% de datos faltantes.
Ahora tenemos que remover dichas estaciones tanto de las series temporales como de los metadatos:
Paso 3: Una vez que tenemos las series temporales y los metadatos de las estaciones que utilizaremos, vamos a cargar los productos de precipitación. Para esto, guardemos las rutas de los productos en dos objetos:
Ahora listemos los archivos que se encuentran en ambas carpetas.
chirps <- list.files(chirps_path, full.names = TRUE)
mswep <- list.files(mswep_path, full.names = TRUE)
head(chirps)
## [1] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900101.tif"
## [2] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900102.tif"
## [3] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900103.tif"
## [4] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900104.tif"
## [5] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900105.tif"
## [6] "../Data/L6_Evaluation/Daily_P_Data/CHIRPSv2//chirps_daily_19900106.tif"
head(mswep)
## [1] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-01.tif"
## [2] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-02.tif"
## [3] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-03.tif"
## [4] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-04.tif"
## [5] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-05.tif"
## [6] "../Data/L6_Evaluation/Daily_P_Data/MSWEPv2.8//MSWEPv280_1990-01-06.tif"Finalmente, podemos utilizar la función rast del paquete terra para importar cada conjunto de archivos:
chirps <- rast(chirps)
mswep <- rast(mswep)
print(chirps)
## class : SpatRaster
## dimensions : 78, 86, 1096 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : 11.625, 33.125, -34.875, -15.375 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : chirps_daily_19900101.tif
## chirps_daily_19900102.tif
## chirps_daily_19900103.tif
## ... and 1093 more source(s)
## names : chirp~00101, chirp~00102, chirp~00103, chirp~00104, chirp~00105, chirp~00106, ...
## min values : 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.0000, ...
## max values : 56.70765, 29.22468, 44.08159, 39.65398, 38.00595, 37.1624, ...
print(mswep)
## class : SpatRaster
## dimensions : 78, 86, 1096 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : 11.625, 33.125, -34.875, -15.375 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : MSWEPv280_1990-01-01.tif
## MSWEPv280_1990-01-02.tif
## MSWEPv280_1990-01-03.tif
## ... and 1093 more source(s)
## names : MSWEP~01-01, MSWEP~01-02, MSWEP~01-03, MSWEP~01-04, MSWEP~01-05, MSWEP~01-06, ...
## min values : 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, ...
## max values : 53.84358, 23.13651, 28.26126, 45.11976, 37.28848, 48.36964, ...Podemos graficar la primera capa de cada set de datos:
Podemos graficar la primera capa de cada set de datos:
Paso 4: Convertir el archivo de metadatos en un objeto espacial para extraer la información de los productos en función de la ubicación de las estaciones.
head(metadata)
## ID LATITUDE LONGITUDE ELEMENT YEAR_START YEAR_END
## 1 BC000068026 -18.367 21.850 PRCP 1932 1992
## 2 BC000068032 -19.983 23.417 PRCP 1921 2022
## 3 BC000068226 -24.017 21.883 PRCP 1958 1992
## 4 BC000068328 -26.050 22.450 PRCP 1939 1992
## 5 BC004612080 -26.467 20.617 PRCP 1960 2022
## 8 BC005847300 -24.667 25.917 PRCP 1925 1997
vct <- vect(metadata, geom=c("LONGITUDE", "LATITUDE"))
print(vct)
## class : SpatVector
## geometry : points
## dimensions : 792, 4 (geometries, attributes)
## extent : 15.05, 32.42, -34.833, -17.817 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : ID ELEMENT YEAR_START YEAR_END
## type : <chr> <chr> <int> <int>
## values : BC000068026 PRCP 1932 1992
## BC000068032 PRCP 1921 2022
## BC000068226 PRCP 1958 1992Podemos usar el paquete maps para agregar el contorno de los paises:
Paso 5: Extraeremos la información correspondiente a la ubicación de las estaciones utilizando la función extract del paquete terra. Esta función nos permitirá extraer los valores de los productos de acuerdo a la ubicación de las estaciones espaciales previamente convertidas en un objeto espacial.
chirps_ts[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## ID 1.0000000 2.000000 3.0000000000 4 5
## chirps_daily_19900101 0.5458916 1.354978 0.0005568612 0 0
## chirps_daily_19900102 0.4041451 0.000000 0.0000000000 0 0
## chirps_daily_19900103 8.3321171 2.423540 0.0000000000 0 0
## chirps_daily_19900104 1.8144772 0.000000 0.0000000000 0 0
mswep_ts[1:5, 1:5]
## [,1] [,2] [,3] [,4] [,5]
## ID 1.0000000 2.000000000 3.000000000 4.0000000 5.00000000
## MSWEPv280_1990-01-01 3.6173794 0.133561790 0.000000000 0.0000000 0.00000000
## MSWEPv280_1990-01-02 0.4801421 0.008661872 0.001714454 0.0000000 0.00000000
## MSWEPv280_1990-01-03 0.8579947 0.005163515 0.000000000 0.0000000 0.00000000
## MSWEPv280_1990-01-04 29.3101654 0.411440879 0.000000000 0.1994748 0.02766745Removemos la primera fila de resultados en ambos data frames.
Podemos agregar el nombre de las estaciones a ambos data frames y convertirlos a objetos zoo.
Es importante destacar que la posición de las estaciones en los metadatos coincide con el orden de las extracciones espaciales. Es decir, la estación en la primera fila de los metadatos corresponde a la primera columna de las extracciones. Es fundamental asegurarse de que el orden de los metadatos sea consistente con el orden de las series temporales de las estaciones. En caso contrario, será necesario ordenarlos de manera correspondiente para garantizar la correcta correspondencia entre los metadatos y las extracciones.
Comparemos el orden de las estaciones entre metadatos y observaciones:
Para evaluar el desempeño de las estaciones, podemos utilizar la función KGE del paquete hydroGOF, como se muestra en el siguiente ejemplo:
Paso 5: En seguida, debemos de desarrollar un loop para evaluar cada uno de los productos con el KGE y sus componentes:
# Vectores vacios para almacenar los resultados
kges_chirps <- r_chirps <- beta_chirps <- gamma_chirps <- c()
for(i in 1:ncol(chirps_ts)){
# Extraer los valores del producto y observados
x <- observations[,i]
y <- chirps_ts[,i]
# Calcular el KGE' y sus componentes
kge_all <- hydroGOF::KGE(sim = y, obs = x, method = "2012", out.type = "full")
kges_chirps[i] <- kge_all$KGE.value
r_chirps[i] <- kge_all$KGE.elements[1]
beta_chirps[i] <- kge_all$KGE.elements[2]
gamma_chirps[i] <- kge_all$KGE.elements[3]
}
# Generar un data frame con los resultados
res_chirps <- data.frame(ID = names(observations), long = metadata$LONGITUDE, lat = metadata$LATITUDE,
KGE = kges_chirps, r = r_chirps, beta = beta_chirps, gamma = gamma_chirps)Podemos visualizar los primeros 10 elementos de los resultados:
head(res_chirps, 10)
## ID long lat KGE r beta gamma
## 1 BC000068026 21.850 -18.367 0.3300503 0.4103995 0.9721374 0.6830971
## 2 BC000068032 23.417 -19.983 0.2477595 0.3962453 0.8472015 0.5781011
## 3 BC000068226 21.883 -24.017 0.3689620 0.4676964 0.9213105 0.6703491
## 4 BC000068328 22.450 -26.050 0.3453056 0.4882832 1.0760378 0.5987656
## 5 BC004612080 20.617 -26.467 0.2843022 0.4233796 1.1143524 0.5917650
## 6 BC005847300 25.917 -24.667 0.2703927 0.3324382 0.7880933 0.7955895
## 7 BC007137570 25.430 -23.120 0.3646196 0.4564381 0.7583529 0.7767167
## 8 BC008382520 21.650 -21.700 0.4369397 0.5456965 1.0085835 0.6674770
## 9 BC008948490 27.500 -21.217 0.4043519 0.5007794 0.8767525 0.6993582
## 10 BC012192590 25.150 -17.817 0.3782020 0.5599314 0.8304033 0.5947725Ahora realizamos lo mismo para MSWEPv2.8.
# Vectores vacios para almacenar los resultados
kges_mswep <- r_mswep <- beta_mswep <- gamma_mswep <- c()
for(i in 1:ncol(mswep_ts)){
# Extraer los valores del producto y observados
x <- observations[,i]
y <- mswep_ts[,i]
# Calcular el KGE' y sus componentes
kge_all <- hydroGOF::KGE(sim = y, obs = x, method = "2012", out.type = "full")
kges_mswep[i] <- kge_all$KGE.value
r_mswep[i] <- kge_all$KGE.elements[1]
beta_mswep[i] <- kge_all$KGE.elements[2]
gamma_mswep[i] <- kge_all$KGE.elements[3]
}
# Generar un data frame con los resultados
res_mswep <- data.frame(ID = names(observations), long = metadata$LONGITUDE, lat = metadata$LATITUDE,
KGE = kges_mswep, r = r_mswep, beta = beta_mswep, gamma = gamma_mswep)Podemos visualizar los primeros 10 elementos de los resultados:
head(res_mswep, warning=FALSE, 10)
## ID long lat KGE r beta gamma
## 1 BC000068026 21.850 -18.367 0.7732305 0.8600063 1.0388965 0.8258930
## 2 BC000068032 23.417 -19.983 0.3683645 0.4052135 1.0246169 0.7888451
## 3 BC000068226 21.883 -24.017 0.7068569 0.7975389 1.0938342 0.8099013
## 4 BC000068328 22.450 -26.050 0.5019049 0.6471822 1.2117271 0.7193046
## 5 BC004612080 20.617 -26.467 0.4202343 0.5808006 1.2333419 0.6744977
## 6 BC005847300 25.917 -24.667 0.5994549 0.7160007 1.0794620 0.7289529
## 7 BC007137570 25.430 -23.120 0.7497016 0.9218014 0.7961388 0.8776282
## 8 BC008382520 21.650 -21.700 0.5602473 0.6129723 0.9501728 0.7972458
## 9 BC008948490 27.500 -21.217 0.8897841 0.9970764 0.9166046 0.9279986
## 10 BC012192590 25.150 -17.817 0.8900011 0.9704940 0.9561920 0.9035116Finalmente podemos comparar el KGE’ y sus componentes:
kges <- data.frame(CHIRPSv2 = res_chirps$KGE, MSWEPv28 = res_mswep$KGE)
boxplot(kges, ylim = c(0, 1),
main = "Desempeño de acuerdo al KGE'", col = c('cyan', 'purple'), ylab = "KGE'")r <- data.frame(CHIRPSv2 = res_chirps$r, MSWEPv28 = res_mswep$r)
boxplot(kges, ylim = c(0, 1),
main = "Desempeño de acuerdo al r", col = c('cyan', 'purple'), ylab = "r")beta <- data.frame(CHIRPSv2 = res_chirps$beta, MSWEPv28 = res_mswep$beta)
boxplot(kges, ylim = c(-1, 1),
main = "Desempeño de acuerdo al Beta", col = c('cyan', 'purple'), ylab = expression(beta))gamma <- data.frame(CHIRPSv2 = res_chirps$gamma, MSWEPv28 = res_mswep$gamma)
boxplot(kges, ylim = c(-1, 1),
main = "Desempeño de acuerdo al Gamma", col = c('cyan', 'purple'), ylab = expression(gamma))